#include <algorithm>
#include <iostream>
#include "3pdm_optimal_irreps.h"
#include "npdm_typedefs.h"

//----------------------------------------------------------------------------------------------------------------------------------------------------------
// Constructor for 3PDM_index object.

ThreePdmBlueprint::ThreePdmBlueprint (int kmax) {

  std::vector< std::vector< std::tuple<int,CRE_DES,ALPHA_BETA> > > ijk_vec, lmn_vec;
  kmax_ = kmax;

  // Build creation 3-indices ordered by site
  get3opsPerms(ijk_vec,CREATION);
  // Build destruction 3-indices ordered by site
  get3opsPerms(lmn_vec,DESTRUCTION);

  printAnalyticalResults();

  buildIrrepOpStrings( ijk_vec, lmn_vec );
  buildSpinDMRGmap();

}

//----------------------------------------------------------------------------------------------------------------------------------------------------------

void ThreePdmBlueprint::printAnalyticalResults()
{
   // Number of sites or spatial orbitals
   int n = kmax_;

   int n2 = n*n;
   int n3 = n2*n;
   int n4 = n3*n;
   int n5 = n4*n;
   int n6 = n5*n;

//   int total = 5*n6/18 - 2*n5/3 + 11*n4/18 + n3/3 - 8*n2/9 + n/3;
   int total = (n-1)*n*(n+1)*(5*n3 - 12*n2 + 16*n -6)/18;
   std::cout << "Irreducible spin-strings = " << total << std::endl;

}

//----------------------------------------------------------------------------------------------------------------------------------------------------------
// Build ordered non-zero permutations of 3 creation/des ops.
// Sorted by site-index

void ThreePdmBlueprint::get3opsPerms(std::vector< std::vector< std::tuple<int,CRE_DES,ALPHA_BETA> > > & pqr_vec, CRE_DES cre_des) {

  std::vector< std::tuple<int,CRE_DES,ALPHA_BETA> > ops(3, std::make_tuple(0,CREATION,ALPHA) );
  std::vector< std::tuple<int,CRE_DES,ALPHA_BETA> > tmp;

  ALPHA_BETA ab;
  for (int i=-kmax_; i<=kmax_ ; i++) {
    if (i==0) continue;
    if ( abs(i)/i == 1) ab = ALPHA;
    else ab = BETA;
    ops[0] = std::make_tuple(abs(i), cre_des, ab);

    for (int j=i+1; j<=kmax_ ; j++) {
      if (j==0) continue;
      if ( abs(j)/j == 1) ab = ALPHA;
      else ab = BETA;

      ops[1] = std::make_tuple(abs(j), cre_des, ab);
      for (int k=j+1; k<=kmax_ ; k++) {
        if (k==0) continue;
        if ( abs(k)/k == 1) ab = ALPHA;
        else ab = BETA;
        ops[2] = std::make_tuple(abs(k), cre_des, ab);
        tmp = ops;
        std::sort( tmp.begin(), tmp.end() );
        pqr_vec.push_back( tmp );

      };
    };
  };

}

//----------------------------------------------------------------------------------------------------------------------------------------------------------

